home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 October / macformat-005.iso / Shareware City / Developers / MacVogl-alpha1PPC / patches.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-07-11  |  9.2 KB  |  490 lines  |  [TEXT/????]
  1. #include "vogl.h"
  2.  
  3. static int    ntsegs, tsegs = 10, nusegs, usegs = 10, 
  4.         ntcurves = 10, nucurves = 10, 
  5.         ntiter = 1, nuiter = 1,
  6.         makeprec_called = 0;
  7.  
  8. static Matrix    et, eu,
  9.         ones =    {
  10.                 {1.0, 1.0, 1.0, 1.0},
  11.                 {1.0, 1.0, 1.0, 1.0},
  12.                 {1.0, 1.0, 1.0, 1.0},
  13.                 {1.0, 1.0, 1.0, 1.0}
  14.             };
  15.             
  16.  
  17. void        premulttensor(), multtensor(), 
  18.         copytensor(), copytensortrans(),
  19.         extract(), iterate(), makeprec(),
  20.         drpatch(), drcurve(), rpatch(),
  21.         replace(), transformtensor();
  22.  
  23. double        addemup();
  24.  
  25. /*
  26.  * defbasis
  27.  *
  28.  *    Specify a basis matrix for a patch or curve.
  29.  */
  30. void
  31. defbasis(id, mat)
  32.     short      id;
  33.         Matrix    mat;
  34. {
  35.         if(!vdevice.initialised)
  36.                 verror("defbasis: vogl not initialised");
  37.  
  38.         copymatrix(vdevice.bases[id], mat);
  39. }
  40.  
  41. /*
  42.  * patchbasis
  43.  *
  44.  *    Specify the two basis matrices for a patch
  45.  */
  46. void
  47. patchbasis(tb, ub)
  48.         long    tb, ub;
  49. {
  50.         if(!vdevice.initialised)
  51.                 verror("patchbasis: vogl not initialised");
  52.  
  53.         copymatrix(vdevice.tbasis, vdevice.bases[tb]);
  54.     copytranspose(vdevice.ubasis, vdevice.bases[ub]);
  55. }
  56.  
  57. /*
  58.  * patchcurves
  59.  *
  60.  *    Specify the number of curves to be drawn in each direction 
  61.  *    on a patch.
  62.  */
  63. void
  64. patchcurves(nt, nu)
  65.         long    nt, nu;
  66. {
  67.  
  68.         if (!vdevice.initialised)
  69.                 verror("patchcurves: vogl not initialised");
  70.  
  71.         if(nt > 0 && nu > 0) {
  72.                 ntcurves = nt;
  73.                 nucurves = nu;
  74.         } else
  75.                 verror("patchcurves: number of patch curves <= 0");
  76.  
  77.     /*
  78.      * Set up the difference matrices
  79.      */
  80.     makeprec();
  81. }
  82.  
  83. /*
  84.  * patchprecision
  85.  *
  86.  *    Specify the lower limit on the number of line segments used
  87.  * to draw a curve as part of a patch. The actual number used varies
  88.  * with the number of curve segments in the "other" direction.
  89.  */
  90. void
  91. patchprecision(tseg, useg)
  92.         long     tseg, useg;
  93. {
  94.         if (!vdevice.initialised)
  95.                 verror("patchprecision: vogl not initialised");
  96.  
  97.         if(tseg > 0 && useg > 0) {
  98.                 tsegs = tseg;
  99.                 usegs = useg;
  100.         } else
  101.                 verror("patchprecision: number of segments <= 0");
  102.     /*
  103.      * Set up the difference matrices
  104.      */
  105.     makeprec();
  106. }
  107.  
  108. /*
  109.  * makeprec
  110.  *
  111.  *    Makes up the two precision matrices for a patch
  112.  */
  113. static    void
  114. makeprec()
  115. {
  116.     double    n2, n3;
  117.  
  118.     /*
  119.      * Find ntsegs, nusegs, ntiter, nuiter....
  120.      * ie. the actual number of curve segments of a tcurve,
  121.      *     the actual number of curve segments of a ucurve,
  122.      *     and the number of times to iterate in each direction.
  123.      */
  124.  
  125.     ntsegs = tsegs;
  126.     ntiter = ntsegs / (nucurves - 1);
  127.     if (ntsegs > ntiter * (nucurves - 1)) 
  128.         ntsegs = (++ntiter) * (nucurves - 1);
  129.  
  130.     nusegs = usegs;
  131.     nuiter = nusegs / (ntcurves - 1);
  132.     if (nusegs > nuiter * (ntcurves - 1)) 
  133.         nusegs = (++nuiter) * (ntcurves - 1);
  134.  
  135.     /*
  136.      * Doing the t precision matrix.....
  137.      */
  138.     identmatrix(et);
  139.     n2 = (double)(ntsegs * ntsegs);
  140.     n3 = (double)(ntsegs * n2);
  141.  
  142.     et[0][0] = et[2][2] = et[3][3] = 0.0;
  143.     et[1][0] = 1.0 / n3;
  144.     et[1][1] = 1.0 / n2;
  145.     et[2][0] = et[3][0] = 6.0 / n3;
  146.     et[2][1] = 2.0 / n2;
  147.     et[1][2] = 1.0 / (double)ntsegs;
  148.     et[0][3] = 1.0;
  149.  
  150.     /*
  151.      * Make the Transpose of eu
  152.      */
  153.     identmatrix(eu);
  154.     n2 = (double)(nusegs * nusegs);
  155.     n3 = (double)(nusegs * n2);
  156.  
  157.     eu[0][0] = eu[2][2] = eu[3][3] = 0.0;
  158.     eu[0][1] = 1.0 / n3;
  159.     eu[1][1] = 1.0 / n2;
  160.     eu[0][2] = eu[0][3] = 6.0 / n3;
  161.     eu[1][2] = 2.0 / n2;
  162.     eu[2][1] = 1.0 / (double)nusegs;
  163.     eu[3][0] = 1.0;
  164.  
  165.     makeprec_called = 1;
  166. }
  167.  
  168. /*
  169.  * patch
  170.  *
  171.  *    Draws a bicubic patch. (ie. all the w coords a 1 and the
  172.  *    basis matrices don't change that)
  173.  */
  174. void
  175. patch(geomx, geomy, geomz)
  176.     Matrix    geomx, geomy, geomz;
  177. {
  178.     rpatch(geomx, geomy, geomz, ones);
  179. }
  180.  
  181. /*
  182.  * rpatch
  183.  *
  184.  *    Draws rational bicubic patches.
  185.  *
  186.  *    Reference: J. H. Clark, Parametric Curves, Surfaces and volumes in 
  187.  *    computer graphics and computer aided Geometric Design.
  188.  *    Technical report No. 221, Nov 1981.
  189.  *    Computer Systems Lab. Dept's of Elecrical Eng. and Computer Science,
  190.  *    Standford University, Standford, California 94305.
  191.  */
  192. void
  193. rpatch(geomx, geomy, geomz, geomw)
  194.         Matrix  geomx, geomy, geomz, geomw;
  195. {
  196.  
  197.     Tensor    S, R;
  198.         Matrix    tmp, tmp2;
  199.     double    xlast, ylast, zlast;
  200.         int    i, j;
  201.     Token    *tok;
  202.  
  203.         if (!vdevice.initialised)
  204.                 verror("patch: vogl not initialised");
  205.  
  206.     /* 
  207.      *  Form S = et . tbasis . Gtensor . ubasisT . euT
  208.      */
  209.  
  210.     if (!makeprec_called)
  211.         makeprec();
  212.  
  213.     mult4x4(tmp, et, vdevice.tbasis);
  214.     mult4x4(tmp2, vdevice.ubasis, eu);
  215.  
  216.     /*
  217.      * Load the geometry matrices into S.
  218.      */
  219.     for (i = 0; i < 4; i++)
  220.         for (j = 0; j < 4; j++) {
  221.             S[0][i][j] = geomx[i][j];
  222.             S[1][i][j] = geomy[i][j];
  223.             S[2][i][j] = geomz[i][j];
  224.             S[3][i][j] = geomw[i][j];
  225.         }
  226.  
  227.     premulttensor(R, vdevice.tbasis, S);
  228.     multtensor(S, vdevice.ubasis, R);
  229.  
  230.     /*
  231.      * Find the last point on the curve.
  232.      */
  233.     xlast = addemup(S[0]);
  234.     ylast = addemup(S[1]);
  235.     zlast = addemup(S[2]);
  236.  
  237.     /*
  238.       * Multiply the precision matrices in.
  239.      */
  240.     premulttensor(R, et, S);
  241.     multtensor(S, eu, R);
  242.  
  243.     if (vdevice.inobject) {
  244.         tok = newtokens(74);
  245.         tok[0].i = RPATCH;
  246.         tok[1].f = xlast;
  247.         tok[2].f = ylast;
  248.         tok[3].f = zlast;
  249.         tok[4].i = ntcurves;
  250.         tok[5].i = nucurves;
  251.         tok[6].i = ntsegs;
  252.         tok[7].i = nusegs;
  253.         tok[8].i = ntiter;
  254.         tok[9].i = nuiter;
  255.  
  256.         tok += 10;
  257.         for (i = 0; i < 4; i++)
  258.             for (j = 0; j < 4; j++) {
  259.                 (tok++)->f = S[0][i][j];
  260.                 (tok++)->f = S[1][i][j];
  261.                 (tok++)->f = S[2][i][j];
  262.                 (tok++)->f = S[3][i][j];
  263.             }
  264.  
  265.         return;
  266.     }
  267.  
  268.     /*
  269.      * Multiply by the current transformation....
  270.      */
  271.     transformtensor(S, vdevice.transmat->m);
  272.  
  273.     /*
  274.      * Draw the patch....
  275.      */
  276.     drpatch(S, ntcurves, nucurves, ntsegs, nusegs, ntiter, nuiter);
  277.  
  278.     /*
  279.      * Set the current (untransformed) world spot....
  280.      */
  281.     vdevice.cpW[V_X] = xlast;
  282.     vdevice.cpW[V_Y] = ylast;
  283.     vdevice.cpW[V_Z] = zlast;
  284. }
  285.  
  286. /*
  287.  * transformtensor
  288.  *
  289.  *    Transform the tensor S by the matrix m
  290.  */
  291. void
  292. transformtensor(S, m)
  293.     Tensor    S;
  294.     Matrix    m;
  295. {
  296.     Matrix    tmp, tmp2;
  297.     register    int    i;
  298.  
  299.     for (i = 0; i < 4; i++) {
  300.         extract(tmp, S, i);
  301.         mult4x4(tmp2, tmp, m);
  302.         replace(S, tmp2, i);
  303.     }
  304. }
  305.  
  306. /*
  307.  * replace
  308.  *
  309.  *    Does the reverse of extract.
  310.  */
  311. static    void
  312. replace(a, b, k)
  313.     Tensor    a;
  314.     Matrix    b;
  315.     int    k;
  316. {
  317.     int    i, j;
  318.  
  319.     /*
  320.      * Not unwound because it only gets called once per patch.
  321.      */
  322.     for (i = 0; i < 4; i++)
  323.         for (j = 0; j < 4; j++)
  324.             a[j][i][k] = b[i][j];
  325. }
  326.  
  327. /*
  328.  * drpatch
  329.  *
  330.  *    Actually does the work of drawing a patch.
  331.  */
  332. void
  333. drpatch(R, ntcurves, nucurves, ntsegs, nusegs, ntiter, nuiter)
  334.     Tensor    R;
  335.     int    ntcurves, nucurves, ntsegs, nusegs, ntiter, nuiter;
  336. {
  337.     Tensor    S;
  338.     Matrix    tmp;
  339.     int    i;
  340.  
  341.     /*
  342.          *  Copy R transposed into S
  343.          */
  344.     copytensortrans(S, R);
  345.  
  346.     for (i = 0; i < ntcurves; i++) {
  347.         extract(tmp, R, 0);
  348.         drcurve(ntsegs, tmp);
  349.         iterate(R, nuiter);
  350.     }
  351.  
  352.     /*
  353.      * Now using S...
  354.      */
  355.     for (i = 0; i < nucurves; i++) {
  356.         extract(tmp, S, 0);
  357.         drcurve(nusegs, tmp);
  358.         iterate(S, ntiter);
  359.     }
  360. }
  361.  
  362. /*
  363.  * iterate
  364.  *
  365.  *    Iterates the forward difference tensor R
  366.  */
  367. static    void
  368. iterate(R, n)
  369.     register Tensor    R;
  370.     int    n;
  371. {
  372.     register int    it;
  373.  
  374. /*
  375.  * Anyone for an unwound loop or two???
  376.  *
  377.  *    for (it = 0; it < n; it++) {
  378.  *        for (i = 0; i < 4; i++) 
  379.  *            for (j = 0; j < 4; j++)
  380.  *                for (k = 0; k < 3; k++)
  381.  *                    R[i][j][k] += R[i][j][k+1];
  382.  *    }
  383.  */
  384.     for (it = 0; it < n; it++) {
  385.         R[0][0][0] += R[0][0][1];
  386.         R[0][0][1] += R[0][0][2];
  387.         R[0][0][2] += R[0][0][3];
  388.  
  389.         R[0][1][0] += R[0][1][1];
  390.         R[0][1][1] += R[0][1][2];
  391.         R[0][1][2] += R[0][1][3];
  392.  
  393.         R[0][2][0] += R[0][2][1];
  394.         R[0][2][1] += R[0][2][2];
  395.         R[0][2][2] += R[0][2][3];
  396.  
  397.         R[0][3][0] += R[0][3][1];
  398.         R[0][3][1] += R[0][3][2];
  399.         R[0][3][2] += R[0][3][3];
  400.  
  401.         R[1][0][0] += R[1][0][1];
  402.         R[1][0][1] += R[1][0][2];
  403.         R[1][0][2] += R[1][0][3];
  404.  
  405.         R[1][1][0] += R[1][1][1];
  406.         R[1][1][1] += R[1][1][2];
  407.         R[1][1][2] += R[1][1][3];
  408.  
  409.         R[1][2][0] += R[1][2][1];
  410.         R[1][2][1] += R[1][2][2];
  411.         R[1][2][2] += R[1][2][3];
  412.  
  413.         R[1][3][0] += R[1][3][1];
  414.         R[1][3][1] += R[1][3][2];
  415.         R[1][3][2] += R[1][3][3];
  416.  
  417.         R[2][0][0] += R[2][0][1];
  418.         R[2][0][1] += R[2][0][2];
  419.         R[2][0][2] += R[2][0][3];
  420.  
  421.         R[2][1][0] += R[2][1][1];
  422.         R[2][1][1] += R[2][1][2];
  423.         R[2][1][2] += R[2][1][3];
  424.  
  425.         R[2][2][0] += R[2][2][1];
  426.         R[2][2][1] += R[2][2][2];
  427.         R[2][2][2] += R[2][2][3];
  428.  
  429.         R[2][3][0] += R[2][3][1];
  430.         R[2][3][1] += R[2][3][2];
  431.         R[2][3][2] += R[2][3][3];
  432.  
  433.         R[3][0][0] += R[3][0][1];
  434.         R[3][0][1] += R[3][0][2];
  435.         R[3][0][2] += R[3][0][3];
  436.  
  437.         R[3][1][0] += R[3][1][1];
  438.         R[3][1][1] += R[3][1][2];
  439.         R[3][1][2] += R[3][1][3];
  440.  
  441.         R[3][2][0] += R[3][2][1];
  442.         R[3][2][1] += R[3][2][2];
  443.         R[3][2][2] += R[3][2][3];
  444.  
  445.         R[3][3][0] += R[3][3][1];
  446.         R[3][3][1] += R[3][3][2];
  447.         R[3][3][2] += R[3][3][3];
  448.     }
  449. }
  450.  
  451. /*
  452.  * Extract the k'th column of the tensor a into the matrix b.
  453.  */
  454. static    void
  455. extract(b, a, k)
  456.     register Tensor a;
  457.     register Matrix    b;
  458.     register int    k;
  459. {
  460.     b[0][0] = a[0][0][k];
  461.     b[0][1] = a[1][0][k];
  462.     b[0][2] = a[2][0][k];
  463.     b[0][3] = a[3][0][k];
  464.  
  465.     b[1][0] = a[0][1][k];
  466.     b[1][1] = a[1][1][k];
  467.     b[1][2] = a[2][1][k];
  468.     b[1][3] = a[3][1][k];
  469.  
  470.     b[2][0] = a[0][2][k];
  471.     b[2][1] = a[1][2][k];
  472.     b[2][2] = a[2][2][k];
  473.     b[2][3] = a[3][2][k];
  474.  
  475.     b[3][0] = a[0][3][k];
  476.     b[3][1] = a[1][3][k];
  477.     b[3][2] = a[2][3][k];
  478.     b[3][3] = a[3][3][k];
  479. }
  480.  
  481. static    double
  482. addemup(m)
  483.     Matrix    m;
  484. {
  485.     return (m[0][0] + m[0][1] + m[0][2] + m[0][3] +
  486.         m[1][0] + m[1][1] + m[1][2] + m[1][3] +
  487.         m[2][0] + m[2][1] + m[2][2] + m[2][3] +
  488.         m[3][0] + m[3][1] + m[3][2] + m[3][3]);
  489. }
  490.